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Abstract 

We apply dynamic mode decomposition (DMD) and proper orthogonal decomposition (POD) 
methods to flows in highly heterogeneous porous media to extract the dominant coherent structures 
and derive reduced-order models via Galerkin projection. Permeability fields with high contrast 
are considered to investigate the capability of these techniques to capture the main flow features 
and forecast the flow evolution within a certain accuracy. A DMD-based approach shows a better 
predictive capability due to its ability to accurately extract the information relevant to long-time 
dynamics, in particular, the slowly-decaying eigenmodes corresponding to largest eigenvalues. 
Our study enables a better understanding of the strengths and weaknesses of the applicability of 
these techniques for flows in high-contrast porous media. Furthermore, we discuss the robustness 
of DMD- and POD-based reduced-order models with respect to variations in initial conditions, 
permeability fields, and forcing terms. 

Keywords: Model reduction, highly heterogeneous porous media, dynamic mode decomposition, 
proper orthogonal decomposition. 

1. Introduction 

In many relevant porous media engineering applications, media permeability can vary several 
orders of magnitude. For example, in flow through fractured porous media, the conductivity within 
fractures can be several orders of magnitude higher than the conductivity within the matrix. Sim- 
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ilarly, shale barriers with very low permeability can affect flow behavior significantly. Because 
of these large variations in coefficient values and the fact that these features (fractures and shale 
barriers) have small spatial dimensions, the direct numerical simulations of these processes may be 
prohibitively expensive. In particular, the resulting large number of degrees of freedom and asso- 
ciated computational cost could inhibit the capability to perform a sensitivity analysis or conduct 
uncertainty quantification studies which require many functional evaluations. As such, construct- 
ing model reduction techniques that can judiciously select the dominant modes corresponding to 
dominant flow features is important in these applications. This enables the derivation of reduced- 
order models with significantly less degrees of freedom while neglecting irrelevant features of 
the physics in order to remain computationally tractable. In this paper, we discuss global model 
reduction techniques for flows in highly heterogeneous media with high contrast. 

Modeling of flow in a high-contrast subsurface requires capturing its long-term dynamics that 
is due to heterogeneous diffusion in the low conductivity regions. In this paper, our goal is to 
develop model reduction techniques that are suitable for accurately predicting the behavior of 
flows in high-contrast media in long time scales. This is motivated by many important applications 
(see [1] and references therein) where the flow response due to low conductivity regions needs to 
be detected in order to identify the location of these regions and their intensity. The sizes of the 
problems involving these features are very large because shale layers or low conductivity features 
can be thin and long and because of the high contrast one needs many grid blocks to resolve these 
functions. In addition, this small-scale fractures may have significant spatial variations of the 
coefficients. Furthermore, these problems usually need to be solved for several initial conditions, 
system's settings, and forcing inputs. To account for all of these, one needs to consider robust 
reduced-order models that enable reliable simplified simulations. 

Several techniques, such as balanced truncation, proper orthogonal decompositions (POD), and 
dynamic mode decomposition (DMD), have been efficiently used for global model reduction, most 
of which involve projection of the original governing equations onto a set of modes. Proper or- 
thogonal decomposition (POD) ||2l- ll2|l constitutes a common technique for extracting the coherent 
structures from a linear or nonlinear dynamical process. This method is based on processing infor- 
mation from a sequence of snapshots and identifying a low-dimensional set of basis functions that 
represent the most energetic structures. These functions are then used to derive a low-dimensional 
dynamical system that is typically obtained by Galerkin projection Q. Schmid Q has 

recently introduced a model-free decomposition technique, namely dynamic mode decomposition 
(DMD), to accurately extract coherent and dynamically relevant structures. This method enables 
the computation, from empirical data, of the eigenvalues and eigenvectors of a linear model that 
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best represents the underlying dynamics, even if those dynamics are produced by a nonlinear pro- 
cess. This technique has been successfully applied for the analysis of experimental Ill5l - l20ll and 
numerical [14, lzi-24] flow field data and has shown a great capability to capture the relevant 
associated dynamics. 

For flows in high-contrast media, long-term dynamics are controlled by the flow in low con- 
ductivity regions. As such, it is essential to keep the modes that capture the local slow dynamical 
behavior in the low conductivity regions when selecting basis functions to represent the solu- 
tion space and derive a reduced-order model by Galerkin projection. As for POD, the selection 
of modes is based on an energy ranking of the coherent structures. However, the energy may 
not in all circumstances be the appropriate measure to rank the importance of the flow structures 
and especially to detect the slow dynamics. Thus, reduced-order models generated by projection 
onto principal components, such as POD modes, unless selected carefully, may be inaccurate be- 
cause the dominant modes identified from a set of snapshots may not necessarily correspond to 
the dynamically-important ones. Moreover, principal components, such as POD, -based reduced- 
order models are often limited in their applications due to their lack of robustness with varying 
flow parameters, initial conditions, and forcing inputs. 

In this work, we apply POD and DMD approaches to flow in highly heterogeneous porous me- 
dia with high contrast to identify the important physical features and derive reduced-order models 
that accurately capture the long-term dynamics of the system. Different numerical examples of 
flows in porous media characterized by varying permeability fields are considered. These perme- 
ability fields include channels and inclusions of high and low conductivity. These configurations 
lead to different types of behavior. The objective is to investigate the capability of POD and DMD 
to capture the main flow characteristics and forecast the flow dynamics response within a certain 
accuracy while reducing the computational cost. In all cases, the DMD-based approach shows a 
better predictive capability and reproduces the flow field with a more reliable accuracy. We observe 
convergence to small errors as the steady-state solution develops when using DMD modes while 
larger errors are obtained when POD modes are considered in the construction of reduced-order 
models. This is mostly due to the DMD's ability to extract the information relevant to the slow dy- 
namics in the system which control the long-time behavior. We also consider parameter-dependent 
problems to investigate the robustness of the POD and DMD modes with respect to variations in 
the initial conditions, permeability field, and forcing inputs. This analysis is motivated by many 
applications where one needs to solve many forward problems corresponding different permeabil- 
ity fields; for instance, when stochastic descriptions are used, such as when the permeability field 
may be subject to uncertainty or multi-phase flow where the permeability is modulated by large- 
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scale mobility. We show that using basis functions generated from a DMD-based approach which 
are determined from few selected samples, one can make accurate predictions of the dynamical 
behavior of the flow in highly heterogeneous porous media. 

2. Problem Formulation and Numerical Examples 

We consider a time-dependent single-phase flow in porous media governed by the following 
parabolic partial differential equation 



where u is the pressure, Q is a bounded domain, / is a forcing term, k(x) is a positive definite scalar 
function that is a function of the spatial location x. k represents the ratio of the permeability over 
the fluid viscosity which is a highly-heterogeneous field with a high contrast (i.e., large variations 
in the permeability). The high and low permeability regions are typically heterogeneous and flow 
within them can be have a complicated structure. 

We consider f2 =]0 l[x]0 1[, u — on dVt, and / is a heterogeneous spatial field repre- 
senting injection and production rates. A finite element mesh is constructed by decomposing the 
domain Q, into N e [ t triangular elements. The variations of the forcing term / in our numerical 
examples over Q is shown in Figure [Q For each cell, the value of / varies randomly between the 
discrete values of -1, 0, and 1. We study the flow behavior under different permeability fields that 
include high and low conductivity regions. We assume homogeneous Dirichlet boundary condi- 
tions in our numerical examples. This model problem is a representation of flow in a reservoir with 
injection modeled by the forcing term. 

The finite element discretization of Equation (HJ yields a system of ordinary differential equa- 
tions given by 



where U and F are vectors collecting the solution values and forcing at the nodes. Here, A = {a^), 
dij = J ■ V0j, M = (rriij), = f n <pi4>j, where 0j is piecewise linear basis functions 

defined on a triangulation of Vt and V denotes spatial gradient. 

Employing the backward Euler implicit scheme for the time marching process, we obtain 



du 
dt 



V(k(i)Vu) + f(x) on Q 



(1) 



MU + AU = F 



(2) 




(3) 
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Figure 1: Spatial variations of the forcing / over the domain Q. 

where At is the time step and the superscript n refers to the temporal level of the solution. 

We analyze two types of the permeability fields, namely those that contain inclusions and 
channels as shown in Figure |2] For each type, we consider two cases of high and low conductivity 
values inside the domain Vt while the background value is kept equal to one. The mesh resolves 
the high or low conductivity regions. For all cases considered in the subsequent analysis, we 
assume the forcing distribution shown in Figure [TJunless stated differently. To get an insight on the 
dynamical behavior of the flow under the different permeability configurations considered in this 
study, we plot the temporal variations of the normalized relative L 2 error between two successive 
solutions | \u l — u l ~ l \ | 2 / | \u l ~ l \ | 2 in Figure|3] As expected, for the cases where channel or inclusion 
permeabilities have low conductivity, the flow field takes longer time to reach the steady state. The 
slow dynamics can be inferred from the eigenvalues of the matrix AtA^ M shown in Figure 
IU We observe that when the permeability has low conductivity regions, there are many eigenvalues 
that are close to 1. The eigenmodes corresponding to these eigenvalues that are close to the unity 
will control long-term dynamics of the flow. The eigenvectors corresponding to these eigenvalues 
simply represent finite element degrees of freedom with support in the low conductivity regions. 
In fact, the number of these modes that are asymptotically close to one when the low conductivity 
parameter approaches to zero is equal to the number of degrees of freedom within low conductivity 
regions. This can be observed by noting that eigenvectors corresponding to small eigenvalues of 
Rayleigh Quotient ^^rar (that corresponds to eigenvalues which are close to the one) consist of 
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(a) Inclusion - low conductivity 



(b) Inclusion - high conductivity 





(c) Channel - low conductivity 



(d) Channel - high conductivity 



Figure 2: Four different configurations of the permeability field. High and low conductivity values inside the inclusions 
and channels are considered while the background value is kept equal to one. 



functions that vary within low conductivity regions while constant in high conductivity regions. 




Figure 3: Convergence to steady-state configuration. 



2.1. Mode decomposition methods 

In the subsequent analysis, we present numerical results to demonstrate the capability of the 
proper orthogonal decomposition (POD) and the dynamic mode decomposition (DMD) techniques 
to reconstruct the flow field. First, we provide a brief description of POD and DMD and their 
implementation. We then discuss the usefulness and effectiveness of these decomposition methods 
to capture the relevant flow characteristics and derive reduced-order models that will be used to 
predict the flow field under different configurations obtained by varying the initial conditions and 
permeability field. 

2.1.1. Proper Orthogonal Decomposition 

A classical way to compute POD modes is to perform a singular value decomposition (SVD) of 
the algebraic operator that maps the states between different realizations, however, this approach 
may have a limited application, especially when dealing with a large mesh size as in direct nu- 
merical simulations (DNS) for which fine meshes (high resolution) imply high computational cost. 
Alternatively, one could use the method of snapshots [25] which allows for a significant reduction 
of the large data sets. In this method, sets of instantaneous solutions (or snapshots) of the flow 
parameters obtained from the DNS are generated and stored in an M — by — N matrix C where 
M and N denote, respectively, the number of grid points and snapshots. Since M ^> N, we seek 
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Figure 4: Eigenvalues of the dynamic process: (a) inclusion type permeability field and (b) channel type permeability 
field.. 

to compute the singular values of C as well as its right singular vector matrix through an eigen 
analysis of the matrix C*C; that is, 



C* C e R 



NxN 



:C*CV, 



and 



tPOD 



-cv t 

Or 



(4) 



where Oi are referred to singular values and V{ are the eigenvectors of the matrix C* C. 

The selection of POD modes is optimal in the sense that the error between each snapshot and 
its projection on the space spanned by those modes is minimized B26I1 . Besides, the square of the 
singular values represents a measure of the energy content of each POD mode and thus provides 
guidance for the number of modes that should be considered in order to capture the relevant physics 
of the system. 

2.1.2. Dynamic Mode Decomposition 

The basic principles and mathematical background of dynamic mode decomposition (DMD) 
are given below following [Q|22]. Over the last few years, this technique has been widely applied 
on experimental [15-2o|] and numerical [3, 2l]-24] flow field data to identify dominant coherent 
structures and help in understanding the underlying physics. These structures can be used to project 
a large-scale problem onto a low-dimensional subspace to obtain a dynamical system with much 
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22B- 



fewer degrees of freedom. For a more detailed description, the reader is referred to Q14J, Il6l 

The DMD method is based on postprocessing a sequence of snapshots to extract the dynamic 
information. Let a snapshot sequence, separated by a constant time step At, collected in a matrix 



8 



Vf ; that is, 

Vf = {vi,v 2 ,v 3 ,-" ,v N } (5) 

where Vj denotes the i^ 1 solution field and the subscript denotes the first element of the sequence 
while the superscript denotes the last element. In Equation ©, these correspond to 1 and N, 
respectively. The basic idea of DMD is to relate the solution field to the subsequent solution 
field Vi+i through a linear mapping A; that is, 

= Avj. (6) 

This assumption leads to a representation of the solution field as a Krylov sequence 

Vf = {v 1 ,Ayi,A ! V" ^-V} (7) 

The objective is to determine the main characteristics of the dynamical process represented by 
the linear mapping A (even if the solution field involves nonlinear aspects). This is performed 
by computing (or approximating) the eigenvectors and eigenvalues of the matrix A. For a very 
large system, these computations may be numerically intractable. Furthermore, for instance in 
experimental fluid simulations, the exact form of the matrix A is not given. As such, an efficient 
and fast numerical approach that approximates well the slow dynamics is useful. 

Assuming that for sufficiently long sequence, the vector can be represented by a linear 
combination of the previous solution fields; that is, 

N-l 

V N =^ a i y i + r ( g ) 

8=1 

or 

v N = Vf - 1 a + r (9) 

where a = {a±, a 2 , • • • , cln-i}* and r is the residual vector. Combining Equations © and ® and 
rearranging the result, we obtain 

A Vf - 1 = Vf = Vf - 1 S + r e^_! (10) 

where = ( ■ • • 1 j is the (N — 1) unit vector and the matrix S is a companion matrix 
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defined as: 



1 



V 



a L 

(1 2 



1 a N _ 2 
1 a^-i J 



(11) 



The unknown matrix S is determined by minimizing the residual r which is obtained by ex- 
pressing the N th snapshot \ N by a linear combination of {vi, v 2 , v 3 , ■ ■ • , Vat_i} in a least-squares 
sense. The minimization problem to determine S is expressed then as: 



S = min || - Vf - X S 



(12) 



To avoid cumbersome notation, we use the same variable for the minimizer as the variable that is 
being minimized. The solution can be determined either using a QR-decomposition of the snapshot 



matrix Vf- 1 ; that is B220, 



S = R 1 Q* V 



N 



(13) 



where the matrices Q and R are obtained from QR decomposition of the snapshot matrix 1 or 



using the Moore-Penrose pseudoinverse of the matrix 1 to obtain a solution Il27l1 



(Vf" 1 )* Vf" 1 



(Vj 



N-l\* v N 
v 2 



(14) 



Once S is computed, the next step is to evaluate its eigenvalues and eigenvectors collected in the 
diagonal matrix D, the matrix X, respectively. 

Finally, the DMD spectrum Xj is obtained by transforming the eigenvalues of S from the time- 
stepper format to the format more commonly used in stability theory (accomplished through a 
logarithmic mapping), and the dynamic modes are computed 4>f MD by weighing the snapshot 
based by the eigenvectors of S; that is, 



\ = log(D,i)/At 



(15) 



I^DMD 



Vf" 1 X; 



(16) 
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where Xj is the j*" column of the matrix X. 

The above algorithm to compute the dynamic modes based on the companion matrix S may be 



ill-conditioned in practical situation. As such, Schmid [|14h proposed a more robust implementa- 
tion. The algorithm includes the following steps. First, we let V^" 1 = U E W* be the singular 
value decomposition of the data sequence Vf -1 . Substituting the SVD representation U E W* 
into Equation (flOl) and multiplying the result by from the left U* and by W E -1 from the right, we 
obtain the following matrix 

U*AU = U*VfWE^ES (17) 

For configuration sets with dim(vj) 3> N, the method of snapshots, as described in Section 12.1.11 
can be used to avoid the computational burden associated with the singular value decomposition 
of large matrices. The matrix U contains the POD modes of the sequence of snapshots 1 , one 
may conclude that the matrix S is obtained from the projection of the linear operator A, which is 
used to approximate the underlying dynamical process, onto a the POD basis. A key advantage 
of the above implementation is the ability to take into account the rank-deficiency of V^" 1 by 
considering a limited basis U given only by the non-zero singular values of E (or by singular 
values above a threshold that can be determined based the amount of cumulative energy that needs 
to be captured). The dynamic modes are computed from the matrix S as follows: 



h DMD 



Uy„ (18) 



where y^ is the j m eigenvector of S, i.e., S y • = /Uj-y^ and U is the matrix collecting the right 
singular vectors of the snapshot sequence \\ 



N-l 



In a recent paper, Chen et al. (12711 proposed an optimized version of the DMD algorithm where 
they employ a global optimization technique to minimize the residual error at all snapshots instead 
of the error at only the last snapshot. They tested the approach over a variety of fluid problems and 
showed its superiority in capturing the relevant frequencies and reproducing flowfields with small 
projection errors. 

2.2. Numerical Examples 

Pre-processing. We solve numerically Equation © over a time interval of [0 8000At] where At = 
0.001. This time interval is observed to be long enough to reach the steady-state solution for 
all cases considered, as shown in Figure |3] Then, we record the first 25 instantaneous solutions 
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(usually referred as snapshots) and collect them in a snapshot matrix as: 



V 



■N 



{vi,v 2 ,v 3 



(19) 



where N is the number of snapshots and the size of the column vectors V j is denoted by M. In these 
numerical simulations, N is taken equal to 25 and M is equal to 1 169 and 1871 for the inclusion- 
and channel-type permeability fields, respectively. The snapshot selection process is determinant 
for the accuracy of the resulting reduced-order model. The use of 25 snapshots is observed to be 
appropriate to reproduce results with acceptable accuracy as will be shown later. 

The POD and DMD basis vectors are computed by applying the numerical approaches as de- 
scribed in the previous section to data obtained from flow simulations in highly heterogeneous 
porous media. For the sake of comparison, we extract and keep the same number of modes for 
POD and DMD in our simulations. 

The POD technique identifies the most energetic structures. The POD is based on an energy 
ranking of the coherent structures obtained by enforcing the orthogonality of the correlation spatial 
matrix. This energy ranking is given by the singular values of the spatial correlation matrix shown 
in Equation ©. Thus, to gain insight into the contribution of each mode to the total energy of 
the system, we define the cumulative energy as E = Ym=i a * an d the cumulative contribution of 
the first j modes as c, = ( X^ =i CjJ /E. We plot in Figure |5] the variations of the normalized 
cumulative energy content with the number of POD modes. Most of the energy is contained in the 
first few modes. Specifically, the first six modes contain more than 99.9% of the total flow energy. 
As such, the following numerical results from the POD and DMD-based approaches use the first 
six modes. 

Approximate solution. We postprocess the snapshot matrix, as described in section 2. 1 , to compute 
the POD and DMD modes and use these modes to approximate the solution field. As such, we 
assume an expansion in terms of the modes <p\ ; that is, we let 



rn 




(20) 



t=i 



or in a matrix form 



U 



n 



u 



(21) 



where $ 





and k can refer to either POD or DMD. 
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Figure 5: Variations of the normalized cumulative energy content with the number of POD modes. 



L 2 projection error. To assess the capability of the POD and DMD modes in capturing the dynam- 
ics involved in the process and enabling good projection subspaces, we project each snapshot onto 
the POD modes, and compute the following inner product 

(Ufa) = (^(x),u(x,t j )) , or a? = ($*$)~ 1 $*U J ' (22) 

where (F,G) = J n (F G) dVL, and define the relative error as the L 2 -norm of the difference 
between the exact and approximate solutions over the exact one; i.e. 

II E(t) \\ 2 = II u{x t) U {x t) || 2 _ 
II u[x,t) \\ 2 

The L 2 projection error is computed for both DMD and POD- based representations. The num- 
ber of modes kept in the expansion given by Equation (|20t is taken equal to six. The corresponding 
results are presented in Figure |6] A low projection error is obtained in the interval of snapshots 
used to compute the modes. This is expected since POD produces the optimal modal subspace 
that contains the largest amount of energy. When using DMD modes to construct the projection 
subspace, a large error is obtained at the first few time steps and recovers as time evolves to reach 
small values. The L 2 projection error increases with time outside snapshots sequence interval and 
converges to steady-state values. We observe that the steady-state error values obtained when ap- 
proximating the solution field in terms of the POD modes are larger than those obtained when 
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Figure 6: Variations of the L2 projection error with time for different permeability configurations. Results are obtained 
using POD and DMD modes. 

projecting the solution field onto the space spanned by the DMD modes. These observations show 
that DMD performs better than POD in extracting relevant dynamic information and then has better 
predictive capability of the long-term dynamical behavior. However, low projection error indicates 
the ability of modal decomposition techniques to compute good projection subspaces, but it does 
not guarantee that a reduced-order model obtained mainly by Galerkin projection of the original 
governing equation on the subspace spanned by the modes will certainly reproduce the reference 
solution. 

Reduced-order model. To obtain the reduced-order model, we use the solution expansion given by 
Equation (|2Q|) . substitute it into Equation (OQ) and project the result onto the space formed by the 
modes as 

L k ,^-V( K (x) Vu)-/\=0, (24) 
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one obtains a set of m ordinary differential equations that constitute a reduced-order model; that 
is, 

a = -($*M$) _1 $*A$a + ($*M$) _1 $*F (25) 

Thus, the original problem with M degrees of freedom is reduced to a dynamical system with 
m dimensions where m <C M. We solve numerically Equations (1251) and compute the approximate 
solution given by Equation (l20l) . Again, six modes (i.e., m = 6) have been used. We consider the 
permeability field shown in Figure |2(a)| and run flow simulations where we vary the number of 
snapshots and compute the errors between the reference and approximate solutions obtained by 
employing POD- and DMD-based approaches. The results are presented in the Tabled] Clearly, 
the selection of snapshots is an important step when applying POD and DMD for model reduction. 
Using only seven snapshots yields large value for the error while considering more snapshots 
improves significantly the forecasting capability of the reduced-order model and leads to small 
errors. The results show that POD and DMD modes computed from 25 snapshots contain relevant 
information of the flow dynamics and then enable a reduced-order model that predicts the flow 
behavior with good accuracy. Increasing the number of snapshots would decrease more the error 
but the aim of this study is to show the potentiality of DMD and/or POD techniques to detect the 
dominant modes that govern long-term dynamics from a small set of snapshots and to forecast the 
evolution of the flow field an acceptable accuracy. 

Table 1 : Variations of the error at t — 5 with the number of snapshots N. Results are obtained using POD and DMD 
modes. 



N 


7 


10 


15 


20 


25 


30 


E \\pod 
E \\dmd 


51.63 % 
37.6 % 


36.8 % 
7.54 % 


19.59 % 
4.69 % 


10.9 % 
3.78 % 


7.9% 
2.8% 


7% 
2.9% 



The variations of the projection error, as defined in Equation (|23T) . with time for different per- 
meability configurations are depicted in Figure [7] Clearly, the large projection steady-state errors 
(more than 80% for the case of the inclusion permeability structure with low conductivity) ob- 
tained when deriving a reduced-order using POD modes shows that POD is able to reproduce well 
the flow field only within the snapshot interval while it obviously fails to predict it outside that 
interval. On the other hand, small projection errors between the reference and approximate solu- 
tions (except those observed at the first few time steps), as shown in Figure Ul demonstrates the 
capability of the reduce-order model obtained by projection the governing equations onto the space 
spanned by the DMD modes to predict the flow field within an acceptable accuracy (at most, the 
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Figure 7: Variations of the Galerkin projection error with time for different permeability configurations. Results are 
obtained using POD and DMD modes. 

error reaches 8%). Clearly, the eigenvalues and eigenvectors of the low-dimensional subspace of 
DMD modes capture the principal dynamics of the flow. In particular, the DMD -based approach 
allows for eliminating the eigenvectors associated with small asymptotically vanishing eigenvalues 
and detect the slow dynamics. The large errors obtained from DMD-based analysis and observed 
at the first few time steps are mostly due to the fact the DMD modes do not contain the eigenvectors 
that correspond to small eigenvalues and which govern the fast decaying dynamics. 

Next, we consider a time- varying permeability field where the permeability is changed in every 



time instant by a mobility function (c.f., two-phase flow simulations 11281301). At this stage, we 
use simplified mobility functions to demonstrate that dynamic modes obtained from DMD can 
be used to accurately predict flow dynamics for different initial conditions and for time-varying 
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Figure 8: Snapshots of the time-varying permeability field, 
permeability fields. As such, we consider the following coefficient: 

Kt(x,t) = k(x) T(x,r(t)), (26) 



where k is the coefficient of the permeability field shown in Figure [2(a)| T(x, r(t)) is defined as 




r(x,r(t)) = t h , (27) 



and n is a circle of a time- varying radius r(t) and center (0,0) and r(t) = 10 t. Figure [8] shows the 
permeability field at three different instants. In this case, the matrix A is time-dependent and then 
it is evaluated at each time step in Equation ([3]). Similar to the previous analysis, we follow POD- 
and DMD-based approaches and derive reduced-order models to investigate their appropriateness 
for time-varying porous media problems. In Figure |9] we plot the variations of the L 2 projection 
and Galerkin projection errors with time. We observe small L 2 projection errors when using POD 
and DMD modes. However, a large error is reached when projecting the governing equations onto 
the space spanned by POD modes to obtain a reduced-order model. This error keeps growing as 
time evolves. On the other hand, a small error is obtained when using DMD modes. This indicates 
the suitability of the use of DMD modes for model reduction of flows in time-varying and highly 
heterogeneous porous media. 

2.3. Parameter-dependent case 

In this section, we investigate the robustness of model reduction techniques with respect to 
moderate variations in the permeability distribution, the contrast, the initial conditions, and the 
forcing inputs. We first consider a permeability field which is represented by a linear combination 
of five different permeability fields with each containing low-conductivity inclusions. These per- 
meability fields are shown in Figure [TOl where each permeability field contains low-conductivity 
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Figure 9: Variations of the L2 projection (represented by dashed lines) and Galerkin projection (represented by solid 
lines) errors with time: time-varying permeability field. Results are obtained using POD and DMD modes. 

inclusions in different locations. This, for example, corresponds to the case where locations of low 
conductivity regions are not deterministic. The coefficient describing the resulting permeability is 
expressed as 



The resulting permeability field obtained for fi 2 , At 3 , /i 4 , fi 5 } = {1,5, 2, 10, 0.1} is depicted in 
Figure [TT1 We compute the POD and DMD modes for each of the five permeability configurations 
(shown in Figure [TO]) and collect them in a global matrix as 



Then, we derive a general reduced-order model, with 5 x m dimensions, by Galerkin projecting 
the governing Equation © onto the space spanned by the POD and DMD global modes and check 
its capability to predict accurately the case shown in Figure [TTJ The temporal variations of the 
Galerkin projection error are plotted in Figure [12] Different initial conditions are considered and 
similar trends are observed. Unlike POD, DMD predicts the flow field with good accuracy. An er- 
ror of 2% is obtained. This error is comparable to the error between the reference and approximate 
solutions obtained when using the DMD modes computed directly for the permeability field shown 
in FigureCEU These results show the robustness of the global DMD modes for developing reduced- 
order models that can be efficiently used to analyze the sensitivity of the dynamical behavior of 
the flow to moderate variations in the structure of permeability field (in terms of distribution and 
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Figure 10: Five permeability fields each containing inclusions with different contrast. 



contrast). 

To investigate further the suitability of POD and DMD modes to model flows in varying and 
highly heterogeneous porous media, we consider the channel-type permeability (shown in Figure 



2(c)) and multiply its coefficient by a smooth positive spatial function; that is, 



k s (x; e; /) = k(x) x (1 + e + sin(27r/x) sin(27r /?/)). 



(30) 



The obtained permeability field for e = 1 and / = 100 is depicted in Figure [T3] We use POD and 



DMD modes generated for the permeability field shown in Figure |2(c)| and employ the Galerkin 
projection to obtain a reduced-order model which is used to predict the flow field resulting from 
the modified permeability field described by Equation (l30l . In Figure Q31 we plot the temporal 
variations of the projection error obtained from the POD- and DMD- based representations while 
varying the value of e. Large Galerkin projection errors are obtained when using POD modes. 
These errors increase substantially as the value of e increases. This indicates that POD-based 
reduced-order model can be only valid for the original configuration considered when computing 
the modes. On the other hand, the DMD-based model reduction approach seems to be much less 
sensitive to variations in the permeability. In fact, it shows a great capability to predict the flow 



field as can be deduced from the small error values shown in Figure 14(a) (about 8% for different 
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Figure 11: Permeability field (a linear combination of five different permeability fields shown in Figure [TUt. 
values of e). 

Next, we modify moderately the distribution and magnitude of the forcing term as shown in 
Figure [15] We use the DMD and POD modes obtained for the forcing case shown in Figured] to 
derive a reduced-order model and predict the behavior of the flow field subjected to the modified 
forcing input. The variations of the L 2 projection and Galerkin projection errors with time are 
plotted in Figure[l6] As expected, smaller errors are obtained from the L 2 projection in comparison 
to those obtained when using the reduced-order model. The use of DMD modes yields small errors. 
This shows the robustness of DMD-based approach to derive a reliable reduced-order model while 
moderately varying forcing inputs. 

3. Conclusions 

In this work, we applied proper orthogonal decomposition (POD) and dynamic mode decom- 
position (DMD) to flow in highly heterogeneous porous media with high contrast to derive a 
reduced-order model. Different numerical examples of flows in porous media characterized by 
highly varying permeability fields were considered. These permeability fields include channels 
and inclusions of high and low conductivity. The long-time dynamics of these flows are due to 
complex changes within low permeability regions. Through our cases, we investigated the capa- 
bility of POD and DMD to capture the main flow characteristics and predict the flow field within 
a certain accuracy. The DMD-based approach showed better capability to reproduce the flow field 
when compared to the POD-based approach. This is mostly due to the DMD's ability to extract the 
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Figure 12: Variations of the Galerkin projection error with time. Results are obtained using POD and DMD modes. 

dynamic information and particularly the modes that govern the long-time dynamics. We also con- 
sidered parameter-dependent problems to investigate the robustness of the POD and DMD modes 
with respect to variations in the initial conditions, permeability field, and input forcing. We found 
that DMD-based approach provides robust basis functions to make accurate predictions of the 
dynamical behavior of flow in highly heterogeneous porous media. 

References 

[1] Y. Efendiev, J. Galvis, E. Gildin, Local-global multiscale model reduction for flows in highly 
heterogeneous media, Journal of Computational Physics 231 (2012) 8100-81 13. 

[2] H. P. Bakewell, J. L. Lumley, Viscous sublayer and adjacent wall region in turbulent pipe 
flow, The Physics of Fluids 10 (9) (1967) 1880-1889. 

[3] L. Sirovich, Turbulence and the dynamics of coherent structures, Quarterly of Applied Math- 
ematics 45 (1987) 561-590. 

[4] A. E. Deane, I. G. Kevrekidis, G. E. Karniadakis, S. A. Orsag, Low-dimensional models for 
complex geometry flows: Application to grooved channels and circular cylinder, Physics of 
Fluids A 3 (10) (1991) 2337-2354. 

[5] G. Berkooz, P. Holmes, J. L. Lumley, The proper orthogonal decomposition in the analysis 
of turbulent flows, Annual Review of Fluid Mechanics 53 (1993) 321-575. 



21 



0.2 0.4 0.6 0.8 1 

Figure 13: Permeability field obtained for e = 1 and / = 100. 

[6] P. Holmes, J. L. Lumley, G. Berkooz, Turbulence, Coherent Structures, Dynamical Systems 
and Symmetry, Cambridge University Press, Cambridge, UK, 1996. 

[7] I. Akhtar, A. H. Nayfeh, C. J. Ribbens, On the stability and extension of reduced-order 
Galerkin models in incompressible flows: A numerical study of vortex shedding, Theoret- 
ical and Computational Fluid Dynamics 23 (3) (2009) 213-237. 

[8] Z. Wang, I. Akhtar, J. Borggaard, T Iliescu, Two-level discretizations of nonlinear clo- 
sure models for proper orthogonal decomposition, Journal of Computational Physics 230 (1) 
(2011) 126-146. 

[9] Z. Wang, I. Akhtar, J. Borggaard, T Iliescu, Proper Orthogonal Decomposition Closure Mod- 
els for Turbulent Flows: A Numerical Comparison, Computer Methods in Applied Mechanics 
and Engineering 237-240 (2012) 10-26. 

[10] I. Akhtar, Z. Wang, J. Borggaard, T Iliescu, A new closure strategy for proper orthogonal 
decomposition reduced-order models, Journal of Computational and Nonlinear Dynamics 
7 (3) (2012) 034503. 

[1 1] A. Hay, J. Borggaard, I. Akhtar, D. Pelletier, Reduced-order models for parameter dependent 
geometries based on shape sensitivity analysis, Journal of Computational Physics 229 (2010) 
1327-1352. 



22 



3.5 




3.5 



3 



2.5 



LU 



1.5 



e=10 E =1 



e = 0.1 



0.5 




0.2 



0.4 



0.6 



0.8 



0.2 



0.4 



0.6 



0.8 



(a) DMD 



(b) POD 



Figure 14: Variations of the Galerkin projection error with time. Results are obtained using POD and DMD modes. 

[12] A. Hay, I. Akhtar, J. T. Borggaard, On the use of sensitivity analysis in model reduction to 
predict flows for varying inflow conditions, International Journal for Numerical Methods in 
Fluids 68 (1) (2012) 122-134. 

[13] M. Ghommem, I. Akhtar, M. R. Hajj, A low-dimensional tool for predicting force decompo- 
sition coefficients and varying inflow conditions, Progress in Computational Fluid Dynamics, 
in press. 

[14] P. J. Schmid, Dynamic mode decomposition of numerical and experimental data, Journal of 
Fluid Mechanics 656 (2010) 5-28. 

[15] D. Duke, J. Soria, D. Honnery, An error analysis of the dynamic mode decomposition, Ex- 
periments in Fluids 52 (2012) 529-542. 

[16] P. J. Schmid, Application of the dynamic mode decomposition to experimental data, Experi- 
ments in Fluids 50 (201 1) 1 123-1 130. 

[17] C. Pan, D. Yu, J. Wang, Dynamical mode decomposition of gurney flap wake flow, Theoreti- 
cal and Applied Mechanics Letters 1 (201 1) 012002. 

[18] P. J. Schmid, Dynamic mode decomposition of experimental data, in: 8th International Sym- 
posium on Particle Image Velocimetry - PIV09, 2009. 

[19] P. J. Schmid, K. E. Meyer, O. Pust, Dynamic mode decomposition and proper orthogonal 
decomposition of flow in a lid-driven cylindrical cavity, in: 8th International Symposium on 
Particle Image Velocimetry - PIV09, 2009. 



23 




Figure 15: Spatial variations of the forcing / over the domain fi. 



[20] F. Lusseyran, F. Gueniat, J. Basley, C. L.Douay, L. R. Pastur, T. M. Faure, P. Schmid, Flow 
coherent structures and frequency signature: application of the dynamic modes decomposi- 
tion to open cavity flow, Journal of Physics: Conference Series 318 (2011) 042036. 

[21] T. W. Muld, G. Efraimsson, D. S. Henningson, Flow structures around high-speed train ex- 
tracted using proper orthogonal decomposition and dynamic mode decomposition, Comput- 
ers & Structures 57 (2012) 87-97. 

[22] P. J. Schmid, L. Li, M. P. Juniper, O. Pust, Applications of the dynamic mode decomposition, 
Theoretical and Computational Fluid Dynamics 25 (201 1) 249-259. 

[23] A. Seena, H. J. Sung, Dynamic mode decomposition of turbulent cavity flows for self- 
sustained oscillations, International Journal of Heat and Fluid Flow 32 (201 1) 1098-1 1 10. 

[24] Y. Mizuno, D. Duke, C. Atkinson, J. Soria, Investigation of wall-bounded turbulent flow using 
dynamic mode decomposition, Journal of Physics: Conference Series 318 (2011) 042040. 

[25] L. Sirovich, M. Kirby, Low-dimensional procedure for the characterization of human faces, 
Journal of Optical Society of America A 4 (3) (1987) 529-524. 

[26] J. Burkardt, Q. Du, M. D. Gunzburger, H. C. Lee, Reduced order modeling of complex 
systems, in: Biennal Conferences on Numerical Analysis, University of Dundee, Scotland, 
UK., 2002. 



24 




Figure 16: Variations of the L2 projection (represented by dashed lines) and Galerkin projection (represented by solid 
lines) errors with time. Results are obtained using POD and DMD modes. 

[27] K. K. Chen, J. H. Tu, C. W. Rowley, Variants of dynamic mode decompositon: Boundary con- 
ditions, koopman, and fourier analysis, Journal of Nonlinear Science DOI: 10.1007/s00332- 
012-9130-9. 

[28] R. Helmig, Multiphase Flow and Transport Processes in the Subsurface: A Contribution to 
the Modeling of Hydrosystems, 1st Edition, Springer, 1997. 

[29] Y. Efendiev, T. Hou, Multiscale Finite Element Methods: Theory and Applications, Vol. 4 of 
Surveys and Tutorials in the Applied Mathematical Sciences, Springer, New York, 2009. 

[30] Y. Efendiev, V. Ginting, T. Hou, R. Ewing, Accurate multiscale finite element methods for 
two-phase flow simulations, Journal of Computational Physics 220 (2006) 155-174. 



25 



